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Abstract 

Background: The activation/repression of a given gene is typically regulated by multiple transcription factors 
(TFs) that bind at the gene regulatory region and recruit RNA polymerase (RNAP). The interactions between 
the promoter region and TFs and between different TFs specify the dynamic responses of the gene under 
different physiological conditions. 

Results: By choosing specific regulatory interactions with up to three transcription factors, we designed several 
functional motifs, each of which is shown to perform a certain function and can be integrated into larger 
networks. We analyzed three kinds of networks: (i) Motifs derived from incoherent feedforward motifs, which 
behave as "amplitude filters", or "concentration detectors". These motifs respond maximally to input 
transcription factors with concentrations within a certain range. From these motifs homeostatic and pulse 
generating networks are derived, (ii) Tunable network motifs, which can behave as oscillators or switches for low 
and high concentrations of an input transcription factor, respectively, (iii) Transcription factor controlled 
adjustable gates, which switch between AND/OR gate characteristics, depending on the concentration of the 
input transcription factor. 

Conclusions: This study has demonstrated the utility of feedforward networks and the flexibility of specific 
transcriptional binding kinetics in generating new novel behaviors. The flexibility of feedforward networks as 
dynamic units may explain the apparent frequency that such motifs are found in real biological networks. 
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Background 

In the last decade or so there have been significant advances in developing new high throughput 
technologies such as microarray, ChiP-Chip or ChiP-seq [2] to uncover large scale cellular networks. 
Understanding these networks has been facilitated by decomposing these networks into smaller parts, 
so-called motifs, with known structure and function [1]. Two major strategies currently exist for 
accomplishing this task. The first one, so called top-down analysis involves mining databases for recurring 
patterns of interacting genes [1,3,4], or in the case of signaling networks, proteins-protein interactions [5]. 
The emerging motifs are analyzed for function, especially the specific roles they might serve in a particular 
network. One interesting result from this work is the discovery of an enrichment in so-called feedforward 
motifs [6] , Fig. 1 (a) , a three gene network. Alon and colleagues [6 9] subsequently demonstrated the 
functionality of these motifs using both theory and experiment. 

Another approach to understanding biological networks is the bottom- up approach [10]. Here knowledge of 
specific protein-DNA interactions and signaling proteins has allowed researchers to synthetically design 
gene regulatory networks. The design of these networks has resulted in the construction of a wide variety 
of network motifs such as oscillators, switches, and logic gates [?, 11, 12]. One of the goals of synthetic 
biology is to design simple networks that can behave as modules. In recent years there has been a 
resurgence of interest in the quantitative properties of gene regulatory networks [?, 13, 14]. In addition, the 
ability to carry out relatively precise measurements of gene activity using synthetic biology techniques [11] 
and single cell measurements [15, 16] has led to considerable progress in this field. 

The top-down and bottom-up approaches are two ways to generate functional modules. Based on a few 
network motifs motivated by bottom-up approaches, we generated several novel functional network motifs 
by exploring additional transcriptional control mechanisms and combinations of feedforward networks. 

One of the principles of gene regulation is regulated recruitment, where transcription factors bind to the 
promoter regions of a gene, recruiting other transcription factors and RNA Polymerase (RNAP) [17]. 
Furthermore, signaling molecules often activate some of these transcription factors through binding or 
phosphorylation. There are multiple ways by which these interactions take place, such as blocking the 



2 



promoter region, twisting the DNA, Histone modification, and DNA looping, etc. [18] . In addition, the 
combinatorial complexity of multiple transcription factors that bind to the same promoter region will also 
determine the transcription rate [19]. Given all this complexity, a very useful approach to model these 
kinds of interactions in gene regulatory networks is the Shea- Ackers approach [20] . In this work, we used 
this method to construct several different kinds of network motifs. 

The network motifs we will discuss can be broadly classified into three categories. (1) Feedforward related 
motifs: obtained by additional feedbacks on the basic incoherent (type I) feedforward network [1]. (2) 
Tunable motif: This subnetwork can either behave as an oscillator or a bistable switch depending upon the 
concentration of a transcription factor. (3) Adjustable gates: This motif can switch between an AND and 
an OR gate depending upon the concentration of a transcription factor. All three network motifs arise 
from multiple transcription factor regulation at a particular gene. In designing these in silico networks, we 
didn't concern ourselves very much with exact mechanism such as DNA looping etc., but rather focused on 
the architecture of the resulting network. 

In the next section, we briefly review the Shea- Ackers approach, which is used to compute the transcription 
rate of genes. Then, we will describe three different kinds of models, amplitude filters, tunable motifs, and 
adjustable gates. The latter two behave functionally very differently depending upon the concentration of 
an input transcription factor. A discussion and summary can be found in the last section. 

Results and Discussion 

Modeling Gene Regulatory Networks Using Shea-Ackers Method 

A common method for modeling protein-gene interactions is by using Michaelis-Menten or Hill function 
kinetics. However, in the models we will describe in this paper, all kinetics will be based on the 
Shea- Ackers formalism [20,21]. This method estimates the probabilities for different transcription states 
from which an overall rate of transcription is derived. The approach is less empirical and more flexible than 
the standard methods and enables one to easily incorporate repressers and activators into the equation. In 
addition the formalism has a relatively straightforward relationship to the stochastic representations of 
gene expression models. Finally, the free energy terms in the formalism can be easily changed by altering 
the promoter region which enables the models based on Shea- Ackers to be tested experimentally. 

We assume that the occupancy of the binding sites on promoters is governed by equilibrium statistical 
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thermodynamics probabilities [20-22]. The probabihties arise from the binding free energies, the free 
energies of interaction between transcription factors bound to adjacent sites, and the concentrations of the 
participating transcription factors and the RNA polymerase. This assumption holds when binding and 
unbinding to the promoter is rapid [21]. Slow binding has been shown to result in stochasticity due to 

operator occupancy fluctuations [23,24], whereas for small concentrations, the individual birth and death 
of proteins introduce noise into the system [25] . 

We illustrate the Shea- Ackers method with a given gene regulated by two transcription factors, Ti and T2. 
We assume that these three proteins bind at three distinct sites on the promoter, providing nine possible 
combinations of binding. Each combination is associated with a free energy that is proportional to the 
probability of this state and can be represented by the following equation [21]: 

/ = exp(-AG//fcsT) [T,r [^2]"^ [R] (1) 

where AG/ is the free energy difference of the bound and unbound states, and the factors [Ti] and [T2] are 
the concentrations of transcription factors, R is the RNAP concentration and ni and 722 are the number of 
monomers which combine into higher order multimers. ks and T are the Boltzmann constant and absolute 
temperature respectively. Prom Eq (1) the normalized probability of a given state is then given by the ratio 

2 eM-^Gf/kBT)[T,r[T,r[R] 
^ E? eM-SGf/knT) [T,r [^2]"^ [R] 

If we assume that the rate of transcription is proportional to the relative probability when the polymerase 

is bound to the gene, then we can partition the states into the polymerase bound state {Zon) and the 

polymerase unbound state (Zoff) [19]. The probability of gene expression is then: 

{Zon + Zqs) ^ ^ 

If a transcription factor is an activator, its free energy to bind with the RNAP will be very low, and hence 
this interaction will be favored. Whereas for an inhibitor, it would be highly improbable for RNAP to be 
recruited for transcription. Therefore, the interactions between the transcription factors and RNAP 
determine the regulatory rules, which we will explore in the next section for several different kinds of 
dynamical networks. 
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Feedforward Related Motifs 

First wc assume that the dynamics himps together transcription and translation into one process. 
Although explicit modeling and experiment has been shown to give rise to interesting effects such as 
protein bursts [25,26] or oscillations as a result of transciption delay [?], we believe that the main features 
of our models are captured by protein-DNA interactions. 

Incoherent Type I Feedforward Related Networks 

High throughput approaches have uncovered several re-occurring motifs termed the feedforward motif [1,3]. 
Fig. 1(a) shows a common transcription factor for two genes, where the third gene is regulated in a 
feedforward fashion. Fig. 1(b) shows a protein-protein interaction in addition to this scheme. Feedforward 
motifs can lead to two types of dynamics depending on the nature of the two signals that converge on the 
third gene. Fig. 2 illustrates a functional representation of the schematic shown in Fig. 1(a) [6,7]. In this 
representation there are three different kinds of interactions: protein degradation, protein synthesis and 
gene regulation. In Fig. 2 the input transcription factor pi modulates the activity of a target gene (G2) 
directly and indirectly through the gene product p2 of another gene (Gl), which is also a transcription 
factor. The interaction between pi and p2 at the promoter region of G2 and their ability to recruit RNAP 
determines the rate of transcription of Gene G2. This type of architecture has been shown [7- 9] to lead to 
two types of dynamics depending on the nature of the regulation that occurs at the target gene G2. If pi, 
P2 are both activators, the gene circuit acts as a low pass filter, i.e. it is able to filter out transient signals 
and transcribe only when the input signal is long lived [6]. If pi, p2 regulate G2 as an activator and 
repressor respectively then the system can act as a bandpass filter, since the delayed response of pi through 
P2 tends to suppress activity of G2 [7] . This has been suggested to be a general mechanism for speeding up 
response times in transcriptional networks [6,9]. Recently, it has also been argued that the steady state 
characteristics of such incoherent feedforward loops could be very important in establishing spatial stripes 
and pulsed temporal expression profiles of transcription factors involved in developmental processes [27] . 

Simple Feedforward 

We assume that pi activates Gl and G2, p2, the gene product of Gl, is recruited by pi, and the protein 
complex pip2 acts as a repressor of G2 (we assume that p2 cannot bind to G2 by itself). These 
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assumptions lead to the following rate equations, 



d[p2\ 
dt 




bi + b2[pi] 



72 [Ps] 



(4) 



dt 



1 + &3 + h[Pl] + b5[Pl][P2 



where the transcriptional rates, and the "lumped" parameters are derived in the appendix. ai,5i represent 
leaky transcription, which is due to the probability that RNAP can bind to the operator region of the gene 
in the absence of a recruiting transcription factor. At steady state, before Gl is saturated with pi, p2 is 
proportional to pi . The transcription rate for G2 (assuming negligible leaky transcription) can be 
approximated as. 



The transcriptional rate rises in proportion to pi for small pi , and falls in proportion to 1 /pi for large 
values of pi. Hence it reaches a maximum at some intermediate value of pi. 

In Fig. 3, we show the steady state values of p2 and p^ with respect to the input pi. We note that the 
concentration of p3 has a maximum for a given value of pi . 

At low input concentration, pi transcribes G2 and Gl, and hence as its input level increases, p2 tends to 
grow. Recruitment of p2 by pi at G2 makes it possible for the piP2 complex to halt further transcription of 
G2. This module is aptly named an "amplitude filter" (originally called a "band detector", [12,28,29], 
since its output is maximal for a specific range of input. Such biphasic response has also been discussed in 
other systems [30,31]. Recently, Ishihara et. al. [27] discussed the band properties of such networks and 
used this to explain pulsed behavior and patterning in Drosophila developmental processes. 

The effect of dimerization 

In this section we will consider the effect of dimerization of p2 , before it binds to the regulatory region of 
G2. Dimerization has been shown to be important to reduce the effects of stochastic fluctuations [32] in a 
feedforward regulatory scheme. Here we discuss the steady state behavior of a network regulated by pi, 
which recruits a dimer of p2, and this complex is a repressor at G2. The only changes that need to be 
made to Eq. |3]are to include the dimerization equations. These modified equations are. 



T oc 



bi'iPi 



(5) 



1 + b2[pi] + bs'lpif 
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dt 

d[p3] 



d[p2] 




- 7i[p2] - 2(A;di[p2]^ - kd2\pd]) 



(6) 



di + rfaPi 



dt 
dpd 
dt 



1 + 4 + d4Pi + dzPiPd 
kdi\p2f - kd2\pd] - kd3\pd\ 



- 12P3 



where in addition to the formation and dissociation of the dimer complex, the dimer can also degrade. 
With these equations it is easy to compute the steady state values of ps with respect to the input pi. In 
Fig. 4, we plot the steady state value of ps as a function of pi in the upper panel, which shows a much 
steeper fall off of the amplitude filter curve, compared to the simple feedforward case discussed earlier. The 
amplitude filter in this case, has a narrower bandwidth, and hence its filtering capabilities of the input 
transcription factor pi arc much more enhanced. The sharp fall off is due to the quadratic dependence of 
the amount of dimer with respect to the input; as the input pi increases the amount of available dimer pd 
to be recruited by pi at G2 increases, thus increasing the amount of repression at G2. Dimerization of pi 
has the effect of increasing both the width and shifting the peak of the amplitude filter curve. 

Effects of mutations at Gl and G2 

We now discuss the effects of two types of mutations at the binding sites at Gl and G2, for the simple 
feedforward with dimerization of p2- Similar qualitative results can be obtained for the simple feedforward 
and mixed feedforward. In general, a mutation at the binding site tends to change the free energy of 

binding of the transcription factor, generally decreasing the binding affinity. In Fig. 4(a), a mutation at Gl 
reduces the ability of pi to bind to it. Hence larger amounts of pi are required to achieve the same 
transcription rate thereby shifting the amplitude filter peak to the right. In Fig. 4(b), a mutation at G2 
reduces the ability of p2 to be recruited to G2. This leads to a slower fall off, since repression does not take 
place very efficiently. The two types of mutations can be used to engineer the shape of the amplitude filter, 
by changing its bandwidth and peak value [?]. 
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Mixed Feedforward Motifs 

In Fig. 5, wc show a slightly different feedforward scheme, which is a functional interpretation of Fig. lb. 
The gene product p2 is shown as an activator of G2, whereas pi is a repressor of G2. Also indicated in the 
figure is a protein- protein interaction, whereby the input pi binds to P2, and targets it for degradation. 

We assume that pi activates Gl leading to the production of p2; P2 can individually bind to G2 and act as 
an activator. Furthermore pi can be recruited by p2 at the operator region of G2, and together this 
complex acts as a repressor, pi binds to p2 and actively degrades it. One possible mechanism by which this 

can occur is if pi labels p2 with ubiquitin molecules for proteolytic degradation [33] . The above 
assumptions lead to the following rate equations, 

d\p2] Ci+C2[pi] 



dt 1 + C3 + C4 [pi 

d[p3] di + d-z [p2 



-7c[Pi]b2] -7i[P2] (7) 
72^3], 



dt I + di + di\p\\ + d^ [P2] + rfe bi] [P2] 

where the transcriptional rates, and the "lumped" parameters are derived in the appendix. In the equation 
for p2, the extra degradation term is due to the protein-protein interaction between pi and p2- For this 
system of equations the plots shown in Fig. 6 show the behavior of the steady state value of p^, 
demonstrating the amplitude filter effect. 

As pi increases, p2 begins to grow and transcribe G2; however, two key factors prevent G2 from being 
continually transcribed with further increases in pi. The first is that P2 is targeted by p\ for degradation, 
and the second is that P2 binds to G2 and recruits pi, which turns off the transcription. 

In both the simple incoherent Type I feedforward, with/ without dimerization, as well as the mixed 

feedforward, we have shown how it is possible to obtain a basic amplitude filter. The property of such an 
amplitude filter can now be further exploited to generate new types of networks when additional feedbacks 
are added to this basic motif. 



Application 1: Time Ordering 

The amplitude filter or " concentration detector" has been shown to perform temporal processing functions 
such as pulse generation [28]. Pulse- like behavior has been simulated in [27], with a series of cascaded 
feedforward loops, which in fact use the amplitude filter property of the networks to generate pulsatile 
behavior. We use a similar idea where a single transcription factor could serve as an input to several 
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amplitude filters, each of which has a different shape. In particular we assume that for each amplitude 
filter module, the peaks of occur at different input transcription factor {piopt) concentration values. 
Then, as the input concentration crosses Piopt^ the amplitude filters get activated in a sequence, depending 
on how far apart piopt are for different amplitude filters [7]. 

A single input can therefore activate several genes in a sequence. One example could be the sequential 
release of different proteins required in a certain developmental process. 

Fig. 8a plots P3 as a function of the input transcription factor concentration, showing a shifted peak for 

the two amplitude filters with respect to the same input transcription factor. In Fig. 8b, the temporal 
profile for the outputs arc shown assuming that the input concentration is ramped up as a function of time. 
This leads to an ordered protein production in time. 

Kashtan et. al. [34], have explored the consequences of multi-output networks regulated by feedforward 
networks. They show through simulations that in some of these cases feedforward loops with a common 
input can regulate genes in a temporal order. Such temporal order can occur in multi-output feedforward 
loop systems such as in the E. colt flagellar synthesis regulation systems, where proteins need to assembled 

in a timed fashion to make up the flagellar basal-body motor [35]. It has also been found in [36], that the 
logic of the program of gene transcription during differentiation in Bacillus Subtilis sporulation involves a 
series of feedforward loops that generate gene transcription in a pulse like manner. 

Application 2: Homeostatic Networks 

Homeostatic networks are important in several biological systems. One example in which homeostasis has 
been shown to occur through integral feedback control is in the chemotaxis network in E. coli ( [37]). Our 
motivation was to design a network using the amplitude filter which would display homeostasis to input 
perturbations. 

Using the output of the amplitude filter module as an input to itself, one can construct the motif shown in 
Fig. 9(a), where the filter module's output ps participates as a transcription factor for gene Gl, whose 
protein product pi serves as the input to the filter module. The feedback could be positive or negative, 
depending on whether the amplitude filter output p^ is an activator or repressor of Gl. In addition, an 
external input protein A can bind as an activator to gene Gl. Depending on the interaction between the 
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input control A, the feedback ps, and RNAP, we obtain different types of behavior. Consider ps as an 
activator of Gl in Fig. 9(a), we assume that the interaction between A, p^ and RNAP is such that we 
obtain an AND gate at Gl, i.e. Gl expresses only when both A and ps are present. The equation that 
describes the additional variable, pi , is given by 



which is the equation for the input of the amplitude filter, and the equation for p^, the output of the 
amplitude filter, is the same as in Eq 6. The steady state plot for the output ps as a function of the input 
A, is shown in Fig. 10(a). 

Figure 10 shows that the steady state values of ps are constant even for large inputs. As the input A of the 
filter module increases, the output ps decreases (assuming that ps is maximal at the initial value of A). 
This is because the input transcription factor concentration moves away from piopt, at which the maximal 
value of transcription occurs. Therefore, this decreases the transcription of the filter module, and since its 
output feeds back as an activator to Gl, its input level tends to decrease the transcription of Gl. 
Essentially the filter module balances the increase in the input A to Gl by decreasing its output. As seen in 
the lower plot of Fig. 10(a), ps stabilizes to an almost fixed value even though the input control A increases 
in time. This is an example of a homeostatic gene network which fixes its response to input transcription 
factor concentrations. The other case shown in the scheme in Fig. 10(b), where is a repressor for Gl. In 
the regulation at Gl, A as an activator and is a repressor. The negative feedback of into Gl 
suppresses the input to the amplitude filter. If A decreases, the input to the filter module p\ also decreases, 
this reduces the transcription rate (assuming that at steady state the value of is maximal) . This then 
lifts the repression from Gl, and the input to the filter module increases, thereby balancing the effect due 
to the reduced input A, and achieving the same steady state as before. The above two circuits produce a 
fixed amount of output proteins, even though the input might vary by a large amount. In the first case it is 
homeostatic to an increase in the concentration value of the input transcription factor, in the second case it 
is homeostatic with respect to a decrease in the input concentration. 

Tunable Motifs 

In this section we discuss an example of a gene network that exhibits oscillatory dynamics at small values 
of an input transcription factor, and bistability for high input values. This network therefore implements 



dpi 
dt 



Ui + U2[A]\p3] 



(8) 



1 + U3 + U4[A] + U5 [P3] + Ue [A] [ps] 
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two different types of motif functionality, depending upon the concentration of an input transcription 
factor. Recently Voigt et. al. [38] discussed a model of a network in the Bacillus sporulation pathway, 
which was shown to exhibit one of two alternatives, i.e. either a bistable switch or oscillatory behavior, 
depending on the environmental conditions. A synthetically constructed network exhibiting 

multifunctionality in E. coli [39] was demonstrated to be able to flip function from an oscillator to a switch, 
by removing a particular interaction. 

There are now several examples of synthetically designed genetic networks such as switches [40,41] and 
oscillators [11,42] that use the common rules of mutual inhibition and activation to achieve a desired 
functional behavior. The network we will describe is motivated by the work of Gardner et. al. [40], where 
the authors designed a toggle switch. In [43], it was discussed how a toggle switch could be converted into 
a relaxation oscillator by suitably manipulating the basic toggle switch network by adding extra regulation. 
Wc further extend this design by introducing a new type of regulation, which involves an external 
transcription factor whose concentration can flip the system function between an oscillator and a switch. 

Fig. 11 shows two mutually repressing genes, Gl and G2. The repression is assumed to occur through 
tetramer binding of each of their gene products, pi and p2, the gene products of Gl and G2, {pi binds to 
G2, and p2 binds to Gl). Hence, if these were the only interacting genes in the network, the system could 
be in one of two stable states, i.e., Gl is fully expressed, and G2 is silent or vice versa. Gene G3 is 
activated by p2, and its product p3 further activates Gl. The feedback of p3 to Gl has the effect of turning 
the bistable switch into a relaxation oscillator [43]. Assuming that initially G2 is ON, which causes G3 to 
get transcribed and ps to grow, ps then activates gene Gl, leading to the growth of pi. Thus Gl switches 
to a high state, which ultimately shuts down gene G2 due to its repressive effects. This leads to a decrease 
in p3 which subsequently turns Gl OFF, thereby completing one cycle of relaxation oscillations. 

Consider now an extra piece of regulation, an external input A, as an activator of gene Gl. A can bind to 
the promoter region and activate Gl, but can also cooperatively bind with the tetramer of p2, which 
represses Gl. ps can also bind cooperatively with the tetramer of P2. which has the effect of repressing Gl. 
The two activators A and p^, however, are assumed to be mutually exclusive, i.e. both A, and ps cannot 
bind together, but each can individually bind to the DNA. As described in the Appendix, these regulatory 
rules make Gl behave like an OR gate with respect to the inputs A and ps. Prom the above regulatory 
mechanisms, the following equations for the protein dynamics emerge. 
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*1 = mi[^]+m2[p3] .g. 

dt l + m3[A]+m4\p2]'^ + m5[A]\p2]^ + m6\p3]+mr[A]\p3] 

dp2 _ ni 

l + n2+n3[pi] 
dpa _ oi + 02 b2] 

dt l + 03 + 04[p2] 



73P3, 



(10) 



We now consider approximating the transcription rate for Gl for small and large values of the control 
transcription factor, A. For small A, 

J, _ "^2b3] 

~ 1 + m4[p2]^ + melpa] + m7[p2]^[P3] 

Gl is activated by ps, and repressed by p2- As discussed earlier ps toggles the bistable switch formed 

between the gene products of Gl and G2. For large A, 

T dl) 

mi ' mi L^'^J 

which is the transcription rate, one would obtain a toggle switch, between Gl and G2. When A is large, it 
is more likely for A to bind to the DNA than its competitor ps and the system is put into the bistable 
regime. The input A can therefore be used to tune the system into a relaxation oscillator, or a switch. Fig. 
12 shows the bifurcation plot for the steady state values of pi as a function of the input transcription 
factor A. The plot shows a subcritical Hopf bifurcation [44,45] for pi ~ 9, and a saddle-node bifurcation at 
Pi ~ 35. As discussed earlier, the system exhibits oscillations for A <~ 11 and bistability for A >~ 35. In 
the right panel of Fig. 12, the time series plots for pi, p2, P3, are displayed, which show oscillations, for the 
case where A=l. In the lower right panel of Fig. 12 the bistable behavior of pi is displayed (for the input 
A = 50), taking one of two values, depending on the initial conditions. 



Adjustable Gates 

The transcriptional interaction between two transcription factors and RNAP has previously been shown to 
generate several instances of logic such as AND, OR, XOR, etc. [19]. Multiple transcription factors can 
regulate the gate properties of a network through their input concentrations. Alon and colleagues [46] have 
studied the gate properties of the regulation of the lacZYA operon in E. coli; and by using a mathematical 
model, their work shows that the regulation can be made to behave as a fuzzy AND, a pure AND, and an 
OR logic gate. In this section we consider gene regulation by three transcription factors which have the 
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property of switching between two different logical functions, depending on one of the input transcription 
factor concentrations. Further we propose that such a circuit with appropriate feedback can be made to 
exhibit other kinds of functionality such as homeostasis or oscillatory behavior. 

Consider gene Gl transcribed by the interaction of three transcription factors, pi, p2 and ps. p^ is the 
control that determines the output logic. We make the following assumptions: binding of all three 
transcription factors and RNAP to the gene is unlikely; individual binding of p-L,P2,P3 and RNAP is also 
unfavorable; we also assume that transcription is not leaky. With these assumptions, transcription occurs 
due to the binding of the following complexes P1P2P, P1P3P, P2P3P ■ The transcription rate, takes the 
following form, 

J, ^ ripip2 + r2PiP3 + r3P2P3 .^2-) 

1 + Tipi + r^p2 + rePa + r'jpip2 + rgpipa + rgP2P3 + rwPiP2P3 

The above formula can be simplified for two cases, i.e. high and low values of the concentration of the 

control transcription factor, p^. For low values of p^, 

T ^ ™ (13) 

1 + npi + r5P2 + rrpip2 

which implies that the transcription is activated only when both pi and P2 are present, implementing an 
AND gate. For high values of ps, the transcription rate is, 

T = r2'pi+r3'p2 

1 + rs'pi + r<)'p2 + rio'pip2 

where r/ = for i = 2, 3, 8, 9, 10. Here the transcription is activated when either pi or p2 is present, 
implementing an OR gate. The control pz is able to switch from an AND to OR gate. In Fig. 13 we plot 
the transcription rates as a function of the input transcription factors, pi and P2- 

The ability to switch from one kind of logical function to another by varying the control ps opens up the 
possibility to use such a motif in a gene network with other interacting genes. Consider Fig. 14(i), a 
regulatory circuit with Gl having three inputs pi,p2 and the control p^, which emerges from a long 
negative feedback loop from its gene product p^ in the following way: The gene G2 constitutively produces 
protein p4, which forms a protein complex with p^. p4 is also a transcription factor for the gene G3, which 
produces protein ps, and hence this is how a feedback into Gl is achieved. We further assume that p4 
binds to the regulatory region of G3 as a tetramer. The transcription factors pi and p2 are assumed to be 
external to the system, and can be either set to a constant or a time dependent value. In our case we shall 
fix the value of p2, but allow pi to fall to low levels, starting at some fixed value. However, the dynamics of 
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the network is symmetric with respect to pi and P2, and hence we could have equally well chosen to vary 
P2 and keep pi fixed. Fig. 15 shows the logical structure for the adjustable gate network. Note that the 
output from the gate determines its behavior. The diflferential equations for the rates of production of the 
various species are, 

d\P5] ^ n [pi] [p2] + r2 [pi] [Pa] + ^3 [P2] Ips] 

dt l + r4[pi] +r5[p2] +r-6[p3] +r-7[pi]b2] + ''sbilbs] + '"9b2]b3] +nobi]b2]b3] 



- ki \P4] IPb] + k2 [C] - bs] , 



dt 

dt 
d\p3] 
dt 



= ki\p4]\p5\ - k2[C], 

= Co- ki [p4] bs] + k2 [C] - \p4] , 



V2\pi]^ 
l + ViipA]"^ 



-Tpsbs], 



For gene Gl, the input p2 is held constant, but the input pi is made to decay from some initial value. Gene 
G2 produces p4 constitutively, and is sequestered by the output of Gl, i.e. p^, into the complex C. Since 
P4 is a transcription factor for G3, its sequestration away from G2 results in a lower value of the protein ps. 
Hence the control ps is changed, the gate properties of Gl switch between AND/OR. The negative 
feedback arising from Gl onto itself is inhibitory, due to the complex formation. In Fig. 14(ii), the upper 
plot shows the steady states of the value of the control ps as a function of the input pi . A supercritical 
Hopf bifurcation is seen to occur at pi ~ 2.5. Fig. 14(ii) lower right-hand plot shows steady oscillations of 
P3. The oscillations arise due to the gate properties of Gl. Initially, when the system is at steady steady 
state, and the inputs pi and p2 are fixed, p3 has a small value which makes Gl behave as an AND gate, 
since both inputs are present, the output p^ is high. There is thus considerable sequestration of p4 due to 
the complex formation. This implies a reduced production of ps, which is consistent with Gl being in the 
AND state. Now if one of the inputs to the system is removed {e.g by making pi decay), since initially the 
gate is in the AND state, the output P5 decreases. This results in the release of p4, which begins 
increasing the transcription of G3; then ps increases, and the gate Gl switches to the OR state. But in 
this state, Gl can be transcribed by p2, and the output p^ increases. This once again results in 
sequestration of p4, and finally reduces the control p^, and in this way we complete one cycle. The system 
therefore switches back and forth between the two states and hence this leads to oscillations. 
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We now describe another application of an adjustable gate. However in this case, the gate properties are 
reversed, i.e., the gate implements an OR gate for low input control transcription factor concentration, and 
an AND gate for high input control transcription factor concentration. We discuss the nature of the 
regulation and its consequences without simulations, since this case is very similar to the previously 
described model. For this case we assume the following regulatory rules: all three transcription factors and 
RNAP bind to the gene; each of the transcription factors pi,P2 and RNAP bind individually to the gene; 
and the complex pi p2 and RNAP can bind to the gene. These assumptions lead to the following rate law, 

J, ^ S1P1P2 + S2PI + S3P2 + S4P1P2P3 .^gX 

1 + S5P1 + S6P2 + S7P3 + SSP1P2 + S9P1P3 + S1QP2P3 + S11P1P2P3 
By inspection it is clear that for low ps, Gl behaves like an OR gate with respect to pi and p2, whereas for 
high ps, it behaves like an AND gate. If we now consider the same network as described above but 
substitute this motif into Gl, then the system shows an almost homeostatic behavior with respect to a 
change in its input. This is easy to understand since initially, when both inputs are present, the output P5 
must be large. Due to the nature of the feedback, this determines the value of ps to be at a low level. The 
system is therefore initially in the OR state. If now one of the inputs is suddenly decreased, since the 
system is initially in an OR state, the output would continue to be high. There would be a small transient 
due to the sudden change in input, but the system would reach a steady state as before. This then allows 
this network to be fairly unperturbed to changes in one of its input transcription factors. 

What could be the function of such networks with adjustable gates? In the first case (with AND gate 
properties for low ps), a sudden drop in a transcription factor concentration could set up oscillatory 

patterns in the network which would then signal the next program to be carried out by the genetic 
network. In the second case (with the gate properties reversed), it is clearly useful to have a homeostatic 
network which works in such a way so as to counteract any sudden changes in input transcription factor 
concentration levels. 

Although we have not found explicit examples for many of the networks we have discussed, which display 
such complex behavior, we believe that such an endeavor is worth exploring. 

Conclusions 

As a first step towards recognizing and understanding large complicated pathways, we have discussed in 
this work the modular design of several functional network motifs. Each of the networks consists of genes 
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which are regulated by multiple transcription factors. The combinatorial regulation was explored in each 
case, and the networks which emerged were found to have very distinct properties. Our modeling 
procedure used the Shea- Ackers method [21], which allowed us to derived the rates of transcription which 
were then used to explore the network dynamics. The networks could broadly be classified into: networks 

which arc derived from incoherent feedforward motifs; and networks which can change their gating 
properties based upon an external input. 

We first discussed the steady state properties of feedforward networks, which can be used as amplitude 
filters. We found that both the Type I simple and mixed feedforward networks led to a similar design, i.e. 
filtering out the input transcription factor concentration, although both networks worked through difii'erent 
types of interactions. Furthermore, we discussed the effects of dimerization for a simple Type I feedforward, 
which has the effect of narrowing the bandwidth of the ampHtude filter. To study the filter characteristics 
of these networks, we simulated the effects of mutations which would change the protcin-DNA binding 
strengths. These generally have the effect of shifting the amplitude filter curve and modulating its 
bandwidth. Furthermore we described how these motifs can be applied to a biological setting. By having a 
common transcription factor as the input to two amplitude filter modules, but with shifted filter 
characteristics, it was possible to obtain a time ordered response of protein production. Homeostatic 
networks emerged if the output of the amplitude filter was fed back to itself. This network was found to be 
resilient in its output to either increasing or decreasing values of an external input transcription factor, 
depending on whether the feedback was assumed to be positive or negative respectively. 

We next described a tunable motif network, where the regulation at one of the genes made it possible for 
the networks to exhibit bistability or oscillatory behavior if one of the external input transcription factor 
concentration was made to increase/decrease respectively. Finally we discussed gate properties of 
regulation at a gene, which can be made to switch between an AND/OR depending on the one of the 
input transcription factor concentrations. 

Methods 
Simulations 

All simulations were carried out using the Systems Biology Workbench (SBW) tools [47]: the network 
designer, JDesigner, the simulation engine Jarnac [48]. Bifurcation diagrams were computed using SBW 
with an interface to MATLAB [49] , and a bifurcation discovery tool [50] . Bifurcation plots were also 
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computed and cross checked using OscillS an interactive bifurcation software package which is linked to 
AUTO [51], and SBW [47]. In aU our simulations the species concentrations are regarded as dimcnsionless, 
whereas the kinetic constants have dimensions of inverse time, with dimcnsionless Michaelis-Menten 
constants. All models are available as Jarnac scripts (supplied in the supplement) which can be easily 
translated to SBML [52] using the JarnacLite tool that is part of the SBW suite [47]. 
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Figures 

Figure 1 - Feedforward motifs extracted from analyzing 

protein-protein and protein-DNA interaction databases, (a) Simple feedforward, the thin arrows indicate 
transcriptional activity, (b) Mixed feedforward; In addition to the transcriptional activity, the thick arrow 
indicates a protein interaction. 



19 



Figure 2 - Interpretations of the feedforward motifs shown in Fig. 1(b) in terms of a reaction network. 

Pi is the input to the system comprising of genes Gl and G2. Gl produces p2, and G2, p3, pi activates G2, 
P2 represses G2. The hue which ends with a dot represents activation, and the hne which ends with a small 
line perpendicular to it represents repression. 

Figure 3 - Steady state values of ps, and p2 for the simple incoherent Type I feedforward network as a 
function of the input factor pi 

Figure 4 - Steady state value of ps as a function of pi for the simple feedforward with dimerization, 
and for the consequences of mutations at Gl and G2. 

Upper plot: The effects of dimerization in a simple feedforward, as compared to a curve obtained for a 
simple feedforward network. Lower Plots: (a) ps as a function of pi for wild (strong binding of pi to Gl) 
and mutant types (weak binding of pi to Gl). (b) ps as a function of pi for wild (strong binding of pip2 to 
G2) and mutant types (weak binding of pip2 to G2). 

Figure 5 - Time ordered production of proteins can be achieved by feeding in a common input into 
two or more amplitude filter modules that have different maximal responses. 

Figure 6 - Sample figure title 

(a): Steady state values of the amplitude filter outputs as a function of the input A, showing the shifted 
output peaks with respect to the input transcription factor concentration, (b): The amplitude filter 
outputs as a function of time showing the time ordering. 

Figure 7 - Homeostatic networks built out of the amplitude filter motifs. 

(a) Feedback, ps from the output amplitude filter positively regulates gene Gl, (b) Feedback, ps from the 
output amplitude filter negatively regulates gene Gl. 

Figure 8 - Homeostatic networks built out of the amplitude filter motifs. 

(a): Steady state values of P3 as a function of the input A. (b): ps, pi and the input A as a function of time. 
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Figure 9 - Interpretations of the second feedforward motifs shown in Fig. 1 in terms of a reaction 
network. 

Pi degrades p2 via a binding reaction, at G2, p2 activates G2, whereas pi acts as a repressor. 

Figure 10 - Steady State response of a mixed feedforward network as a function of input factor pi. 

Figure 11 - A cartoon of the tunable motif which shows three interacting genes: Gl, G2 and G3, and 
their gene products pi,p2,P3 respectively. 

Regulation at Gl includes positive regulation by an external transcription factor A and p3, and repression 
by p2. The transcription factor A for large concentrations can out compete ps in binding to Gl. The effect 
of this is to break the feedback between the two antagonists which make up the toggle switch, and hence 
the network exhibits bistability for large A. 

Figure 12 - Bifurcation diagram and time-series plots for the tunable motifs. 

The left figure shows the bifurcation plot for the steady state value of pi, as a function of the transcription 
factor concentration A. A subcritical bifurcation occurs at A ~ 9, and oscillations occur for A <= 11. At 
A ~ 35, a saddle-node bifurcation occurs which gives rise to the appearance of two stable steady states. 
Figure (a) shows time series plots for pi, P2, Ps, which show oscillations, for the case where A = 1. Figure 
(b) shows the final steady state of pi reaching either of two values, depending on the initial conditions. 
This is due to the bistable behavior of the circuit. The parameter for this case is A = 50 

Figure 13 - The figures show the transcription rate as a function of the transcription factor 
concentrations, pi,P2> for fixed values of the control transcription factor p3. 

(a), p3 = 1, represents an AND gate, (b), ps = 10, represents an OR gate. 
Figure 14 - An application of adjustable gates 

Panel (i) shows the schematic of the network while Panel (ii) shows the simulations results. In (ii), the 
upper plot shows the bifurcation plot of ps, as a function of the input pi. A supercritical Hopf bifurcation 
is seen to occur at pi ~ 2.5. The lower panel in (ii)(a) shows the decay of pi, while the lower panel in 
(ii)(b) shows the system reaching steady oscillations of ps, after an initial transient. 
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Figure 15 - A logical diagram of the schematic shown in Fig. 14. 

The triangle symbol represents a Boolean inverter. The action of the Boolean gate switches the AND-OR 
gate between AND and OR functions 



A Appendix 1 

In this section we demonstrate the equivalence of the thermodynamic (Shea- Ackers) approach with that of 
the enzyme kinetic approach to derive the equations for the feedforward network. Although this is well 
known we provide this derivation explicitly to show how a particular choice of interaction can be used to 
determine a reaction scheme. The reaction scheme is required if a stochastic simulation is to be performed. 
Prom Fig. 2 we see that pi is an activator for Gl. RNAP can also bind to Gl in the absence of pi, but at a 
lower rate. As described in the main text, such binding leads to "leaky transcription" . At the operator 
sites of G2, pi can bind and recruit RNAP to the operator site of the gene. However, pi also recruits p2, 
the gene product of Gl, and in this state, transcription is repressed. We also assume that p2 binding to 
G2, occurs only as recruitment, and cannot occur by itself. Prom these assumptions, we can draw up the 
following truth tables for the transcriptional regulation at Gl and G2. 
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Table 1: Logic for the transcriptional regulation underlying simple feedforward genes. In the case where 
transcription is unlikely to occur the rate is denoted by x and assigned to zero in the calculations. 

In the above table and the subsequent tables the terms ai,(3i, are related to the free energies of binding 
through exp The logic in Table 1 can be translated to the rate of transcription by computing the 
fractional probability of each of the genes Gl and G2 being bound by RNAP. The transcription rates, 
TrGi,TrG2, are proportional to the probability of occupancy of RNAP, which can be computed from the 
above two tables to be, 

Tr Of ^1 [P]+«3 [Pi] [P] . . 

1 + ai [P] + a2 [pi] + as [pi] [P] ^ ' 
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Table 2: Reaction Scheme for the feedforward network. 

f34 [P]+f}6 [Pl] [P] 



TrG2 oc 



1 + /34 [P] + 132 [Pi] + /Ja [Pi] [P2] + /36 [Pi] [P] 



The transcription rates are then used in the rate laws described in Eq. 4. We now compute the above 
result using the reaction kinetics approach. Assume that for the gene Gl, there are 4 states: Gl-unbound 
or free, Glpi-bound by p\, GlP-bound by RNAP, GipiP-bound by p\ and RNAP. Then it follows that, 

[Gl] + [Glpi] + [GIP] + [GlpiP] = 1. (18) 

For gene G2, there are 5 states: G2-unbound or free, G2pi-bound by pi, G2P-bound by RNAP, 
G2pip2-bound by pi and p2 and G2piP-bound by pi and RNAP, which leads to, 

[G2] + [G2pi] + [G2P] + [GlpiP] + [Glpips] = 1- (19) 

The following reaction scheme defines the network. Prom the above reaction scheme, we obtain at 

thermodynamic equilibrium, 

[Glpi] = fci [pi] [Gl], [GlpiP] = k2p [Glpi] [P\, [GIP] = kip [Gl] [P] (20) 
[G2pi] = fca [G2] [pi], [G2pip2] = ki [G2pi] [pa], [G2P] = k^p [G2] [P], 
[G2piP] = kip [G2pi] [P] 
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where in the above equations we use the equihbrium constants, which are ratios of the forward to backward 
rates, e.g k\ = etc. Using Eqns. 17, 18 & 19, the ratio of genes bound by RNAP, and hence the 

"■16 

transcription rates can be evaluated to be, 

Tr oc fcip P + ki k2pPi P 

''^^ °^ 1 + fci pi + ki k2p piP + kipP ^ ' 

ksp P + k2 kippi P 



l + k2,p P + k2 pi+ ki kz pi p2 + /c2 k4p pi P 



which is functionally the same in form to Eq. 17. Both methods are equivalent as the only physical 
requirement is thermal equilibrium for these reactions. Although the statistical approach is more intuitive, 
the equilibrium approach allows us to define a reaction scheme, which can be described in terms of 
measurable kinetic constants. 



B Appendix 2 

In each of the subsections below, the transcriptional rates for the genes, which are described by the truth 

table for the transcriptional rules, are derived. The parameter values used for the simulations arc also 
provided. The parameters used in the equations are lumped, in the sense that we group together terms 
which depend on the same variable, e.g [pi] [P] + a2 [pi] = vi p\ etc. We also assume for simplicity that 
P = 1. 



B.l Feedforward Networks 

B.1.1 Simple Feedforward 

Table 1, in Appendix I describes the regulation at Gl and G2. The parameter values used in Eq. 4 are. 
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Table 3: Parameters values for Eq. 4. 



B.l. 2 Mixed Feedforward 
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Table 4: Logic for the transcriptional regulation underlying mixed feedforward network genes, Eq. 7. 
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Table 5: Parameters values for Eq. 7. 



B.1.3 Consequence of Dimerization for a Simple Feedforward and the effects of Mutations at Gl, G2. 
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Table 6: Parameters values for Eq. 6. 



The upper plot in Fig 4, for the simple feedforward (monomer) we used the same parameters values as in 
Tables. For the lower two figures, the mutations were carried on the parameters 62, d^. Their values are 
given by, panel A: 62 = 0.001-wild type, 62 = 2 x 10~^-mutant, panel B-.ds = 0.003-wild type, 
(^5 = 5 X 10~^-mutant 

B.1.4 Application 1: Time Ordering 

The amplitude filter parameters used for the simulations in Fig. 6 are the same as in Table 3, except for 
the following changes given by, 



25 



Filter # 
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Table 7: Parameters values for Fig. 6. 



B.1.5 Application 2: Homeostatic Networks 

The regulation at Gl in Fig. 7 (a) is assumed to be, 



A p3 P Rate 

~0 F 

1 pi[A] 

1 P2[P3] 

1 1 P3[A]\p3] 

1 P4[P] 

10 1 X 

011 X 

1 1 1 P7[A]\ps][P] 

Table 8: Logic for the transcriptional regulation underlying Gl for the homeostatic network. 



r,^^ ^ P4[P]+P7[A][ps][P] 

1 + PI [A] + p2 bs] + P3 [A] [P3] + P4 [P] + P7 [A] [P3] [P] 

The parameter values used in Eq. 8, to generate Fig. 8 are given by, and the parameters used for the 
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Table 9: Parameters values for Eq. 7. 



amplitude filter are the same as in Table 4. 



B.2 Tunable Motifs 



Trci oc 
TrG2 oc 
TrG3 oc 



ee[A][P] + sr\p3][P] 



1 + si[A] + e2[P2*] + es[A]\p2^] + s^lps] + S5\p2^]\ps] + ee[A][P] + e7\ps][P] 

C2[P] 



l + C2[P]+Cl[Pl1 

V2[P]+V3\P2][P] 
1 + Vl\P2] + V2[P] + V3\P2][P] 



The parameter values in Eq. 8, are given by, 
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Table 10: Logic for the transcriptional regulation underlying the Tunable Motif network genes, for Gl, G2 
and 03 respectively. 
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Table 11: Parameters values for Eq. 8. 



B.3 Adjustable Gates 
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Table 12: Logic for the transcriptional regulation underlying the adjustable gates, i.e. regulation of Gl, for 
the two cases considered in the main text (Eq. 12 & Eq. 15 respectively). The Table on the left refers to 
the case for which Gl behaves as an AND gate for low pa, and an OR gate for high ps. The table on the 
right is for the opposite case. 



Zonl + Zoffl ' 

Zonl = Os bl] bs] [P] + O9 [P2] ba] [P] + 010 bi] b2] [P] , 
Zoffl = 1 + 6»ibl] + 6*2 b2] + 6*3 bs] + 6'4bl]b3] + 6'5b2]b3] + 6'6bl]b2] + 6'7bl]b2]b3], 
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rp Adjustgate2 ^o«2 
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Zoff2 = l + +^?2b2] +t^3[P3] + t^4bl][P3] +^5b2][P3] +^?6bl][P2] +^7bl][P2]b3]- 



p 


p 


Rate 








1 


1 










1 


MP] 








Mnp] 



[pa]'' P Rate 

~~0 1 

1 xi[p4^ 

1 X2[P] 

X3b4]'[^'] 



Table 13: Logic for the transcriptional regulation for genes G2 & G3 in Fig. 14. 
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Table 14: Parameters values for Eq. 15. 
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